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Abstract 

The problem of calculating the stability of steady state solutions of differential equations is 
treated. Leading eigenvalues (i.e., having maximal real part) of large matrices that arise from 
discretization are to be calculated. An efficient multigrid method for solving these problems is 
presented. The method begins by obtaining an initial approximation for the dominant subspace 
on a coarse level using a damped Jacobi relaxation. This proceeds until enough accuracy for 
the dominant subspace has been obtained. The resulting grid functions are then used as an 
initial approximation for appropriate eigenvalue problems. These problems are being solved 
first on coarse levels, followed by refinement until a desired accuracy for the eigenvalues has 
been achieved. The method employs local relaxation on all levels together with a global change 
on the coarsest level only, which is designed to separate the different eigenfunctions as well as 
to update their corresponding eigenvalues. Coarsening is done using the FAS formulation in a 
non-standard way in which the right hand side of the coarse grid equations involves unknown 
parameters to be solved for on the coarse grid. This in particular leads to a new multigrid 
method for calculating the eigenvalues of symmetric problems. Numerical experiments with a 
model problem demonstrate the effectiveness of the method proposed. Using an FMG algorithm 
a solution to the level of discretization errors is obtained in just a few work units (less than 10), 
where a work unit is the work involved in one Jacobi relaxation on the finest level. 

^Supported in part by the National Aeronautics and Space Administration under NASA 
Contract Nos. NAS1-18107 and NAS1-18605 while the author was in residence at ICASE, NASA 
Langley Research Center, Hampton, VA 23665, and in part by the Air Force Office of Scientific 
Research, United States Air Force under Grant AFOSR-86-0127. 
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1 Introduction 


Stability analysis plays an important part in hydrodynamics [4] and in other physical sciences [6]. 
Basically, one is interested in the growth (decay) rate of infinitesimal perturbations to steady state 
solutions. In rectangular geometries with simple enough steady state solutions one can use Fourier 
transform techniques to reduce the resulting stability problem to an eigenvalue problem for an 
ODE with parameter(s) [4]. General geometries or complicated solution structure do not permit 
this and eigenvalue problems in two and three space dimensions need to be solved. 

A problem related to the stability of a given steady state solution is the problem of transition, 
when a nonlinear problem with a parameter is given and the value of that parameter at which 
steady state solutions loose their stability is required. In fact, in cases of marginal stability this is 
the relevant problem instead of finding the stability of a given steady state solution, since in such 
cases the discretization (in space) might have shifted the eigenvalue from, say, a stable region to 
an unstable one. 

Steady state solutions are often obtained by time marching algorithms. In such cases the 
stability to small perturbations can be studied during the solution process. However, efficient 
multigrid solvers for steady state problems do not involve usually any time marching. The solver 
is designed for solving the steady state problem directly. This allows obtaining solutions which are 
physically unstable. A way to recover the information about the stability of a given steady state 
solution is therefore needed. 

One way of obtaining that information is to use time marching for the time dependent problem 
obtained by linearizing the original problem, around the steady state solution, using a random 
initial condition. The discretization has to be such that it preserves stability. That is, if a certain 
component is decaying oscillating or growing, the same should happen for the discretization. Crank- 
Nickolson time stepping has this property. 

The effectiveness of such a method depends on the closeness of the dominant eigenvalues to 
the imaginary axis. The closer they are the more time steps are needed to determine the growth 
rate of small perturbations to the steady state solution under study. For symmetric problems the 
situation is much better than the non-symmetric ones, since the time steps can be taken as larger 
and larger for larger times. In the non-symmetric case the time stepping depends also on the 
imaginary part of the dominant eigenvalues, and this can be much larger than the real part. In 
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other words, the time stepping can depend on the oscillation rate rather than the decay rate which 
may be much smaller. An improvement of the above method is to use time marching until the 
subspace of dominant eigenvalues has been determined to a desired accuracy, then obtaining from 
it the set of eigenfunctions and eigenvalues as in [7]. 

The above approach is inappropriate when the steady state solution has been obtained by an 
efficient multigrid solver since it is more expensive than getting the steady state solution using a 
multigrid method. Moreover for transition problems it can be very time consuming and a different 
approach should be used. 

i In this paper an alternative way for calculating the stability of steady state solutions by solving 

directly, using a multigrid algorithm, an appropriate eigenvalue problem. This approach could be 
modified for transition problems. The difference is only in a global step (explained later) that is 
performed on the coarsest level. The exact treatment of transition problems by multigrid, together 
with some physical examples will be described elsewhere. 

Since the eigenvalue with maximal real part is required, the multigrid processes must not 
contradict this requirement. The relaxation should be such that it will relatively damp all the non- 
dominant eigenfunctions. Jacobi or a damped Jacobi relaxation achieves this purpose. Kaczmarz 
on the other hand, should not be used since it damps more the eigenfunctions whose eigenvalues are 
away from zero. As in [2] a single grid method which involves a relaxation method and a step for 
updating the eigenvalue is constructed first. In non-symmetric problems the relaxation (Jacobi- ) 
couples the eigenfunctions that correspond to complex conjugate eigenvalues. Therefore, the step 
! for updating the eigenvalue, (the global step) performs also a projection onto one of the subspaces 

of the two conjugate eigenvalues, say, with non-negative imaginary part. 

When the dominant eigenvalue is not separated enough from the others, in the sense that their 
real parts are close, one has to include more eigenfunctions and to solve for the dominant eigenspace 
which is of higher dimension. 

The coarsening strategy for the multigrid method is done using an FAS formulation since the 
eigenvalue problem is a non-linear one. However, since the global change done on the fine level 
affects the coarse grid right hand side significantly in this problem, the standard FAS should be 
modified to include this change of the right hand side. The resulting FAS is such that the right 
hand side of the coarse grid equations involves some parameters to be solved on that level. The 
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resulting algorithm has fast algebraic convergence rates. Moreover, it leads to a new multigrid 
method for symmetric eigenvalue problems. 

The algorithm is performed in an FMG version, where the problem is solved first on coarse levels 
using the approximation obtained there as an initial guess for the solution of finer grid problems. 
On each level a fixed number of basic multigrid cycles are performed, yielding a solution to the 
level of discretization errors in just a few work units (about 10), where a work unit is the work 
involved in one Jacobi relaxation on the finest level. 

Numerical experiments with a two dimensional model problem are given. The results demon- 
strate the effectiveness of the proposed method. Solution to the level of discretization errors is 
obtained with the 1-FMG algorithm (which uses one basic multigrid cycle per refinement and cost 
about 10 WU). 

2 The Problem 

Let A be a linear real valued operator. We wish to compute its eigenvalues with maximal real part. 
That is 


A$t = W with Real(A) maximized. (2.1) 

In general, for asymmetric operators, the eigenvalues and eigenfunctions are complex. The above 
eigenvalue problem can then be written in terms of real quantities as 


(A-\I)<f> = 0, \\<f>\\ 2 = 1 


( 2 . 2 ) 


or as 


( A- xi ni \ ^ 

[ -til A-XI J\rP y 



m 2 = \m 2 = i 


(2.3) 


depending on whether the dominant eigenvalue is real or complex. In this notation A = Real( A), 

H = Imag( A), <f> = Realft), xj) = Imag( S'). 

Since the operator A is real its eigenvalues come in pairs of complex conjugates. Hence, if (A, /z) 
is a solution to (2.3), also (A, — h) is a solution. Moreover, each of these eigenvalues is double. If 
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(cf>,tp) T corresponds to (A ,/z) then so does (tp,—<f>) T . Similarly, (rp,<f>) T and (<j>, — rp) T correspond 
to (A / 

3 Single Grid Algorithm 

One way to get functions in the subspace spanned by the eigenfunctions with eigenvalue of maximal 
real part is to simulate a time dependent process, that is, solving approximately the initial value 
problem 

U t = AU (3.1) 

17(0) = random (3.2) 

See (I. Goldhirsch et al 1987). The dominant part of the solution at large time is a linear combi- 
nation of the desired eigenfunctions. This leads us to the following. 

3.1 Initial Approximation 

The initial approximation for the eigenvalue problem is obtained from the following process applied 
repeatedly, 

(a) Relax using a damped Jacobi relaxation AW = 0. 

(b) W*-W/\\W\\. 

The stopping criterion for the process is as follows. Let </>i = W, <f>j = A4>j-i {j > 1), and let 
be an nxn matrix whose (ij) element is < >. Define 1 to be the minimal positive integer 

for which det — ► 0. The dominant subspace A is of dimension 1, and the process can be 

stopped when det £>(* +1 ) is small enough. The dimension of A determines how many eigenfunctions 
are to be simultaneously solved for in order to obtain fast convergence. 

After an initial approximation for the dominant subspace has been obtained a solution of the 
eigenvalue problem can be started. The algorithm consists of two steps which are performed succes- 
sively. The first is a local relaxation process on equation (2.2) or (2.3) depending on the eigenfunc- 
tion. Its role is to damp the error in the components corresponding to non-dominant eigenvalues. 
The analogy of damped Jacobi relaxation with explicit time stepping suggests that using it with 
an appropriate under-relaxation parameter may serve our purpose. The local relaxation is then 
followed by a step referred to as the global step which is explained in detail next. 
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The Global Step 

(a) dim M = 1. In this simple situation A can be updated according to 

A - < /*’*^ . * - H\m (3-3) 

< <P,9 > 

(b) dim M = 2. Assume for simplicity of exposition that the dominant subspace is spanned by a 
pair of complex conjugate functions. In that case the approximation we get at an intermediate step 
will be a linear combination of the four eigenfunctions mentioned in section 2, and the eigenfunctions 
corresponding to other A’s. The latter converge to zero in the relaxation process (relative to these 
four), while these four corresponding to the maximal A are coupled by the relaxation process. If our 
approximate solution is in the subspace spanned by the two eigenfunctions corresponding to one of 
the /i’s, say, the non-negative one, which belong to maximal A, then A and n can be calculated by 
certain inner products. Since in general this is not the case, we have to perform a projection onto 
the two dimensional subspace that corresponds to one of the /*’ s, say, the nonnegative one. 

Let (4>, ip) T be an approximate solution of (2.3) which lies mainly in the relevant four dimensional 
subspace. A projection into a two dimensional subspace which correspond to (A, fJ.) is to be done: 
Each of the two components of an eigenfunction is approximately a linear combination of the current 
<p and ip. That is, 


, \ 


f , \ 

<p ) 

at ix ai 2 

1 * 



_ 

J 

^ <*21 022 j 

\*J 


(3.4) 


Substituting this in (2.3) and taking inner products of each of the two equations <p and ip, we get 
a 4x4 system given below 




( \ 


C-XD 

t*D 

( * 

= 

, D 

C-XD } 

V* J 

l°J 


(3.5) 


||aii<£ + ai2V’|| 2 = 1 


(3.6) 


||0!21^ + «22V , I| 2 — 1 


(3.7) 
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where C and D are given by 


< A<f>,<f> > < Atf),<f> > 

K < A<f>,ip > < Atp,tp > J 


(3.8) 


D = 


<4>,<t>> <tp,<f>> 
<<j>,rp> <<p,rp> 


(3.9) 


and Oj = (an,ai 2 ) T , a 2 = (o' 2 l> 0! 22 ) T - This is equivalent to the following (complex valued) 
generalized 2x2 eigenvalue problem 


o 

II 

el 

<a. 

1 

O, 

(3.10) 

||i?ea/[a(<£ + »V>)]|| 2 = 1 

(3.11) 

\\Imag[a[<t> + ixl>)]\\ 2 = 1 

(3.12) 


where p = A ± ip. The solution of the above eigenvalue problem gives A = Real(p), p = Imag(p) 
and Oj = Real (a). a 7 = I mag [a). Under our assumption the matrix D is nonsingular since <f> and 
i p are linearly independent. Hence, the above 2x2 eigenvalue problem is a standard one, and has 
two eigenvalues. 

(c) dim "H = n. 

Let M be spanned by the set of real valued functions (<f> i, . . . , <j> n ). A projection into the subspace 
corresponding to non-negative imaginary part . . . , <f> n ) T <— E{<j> j, . . . , <f> n ) T is described. Let C 
be an nxn matrix whose (ij) element is < A<t>i><j>j > , and D is an nxn matrix whose elements are 
< <f>iy<f>j >. D is non-singular since the set of <f>'s is linearly independent. The projection step is 
done by solving the nxn eigenvalue problem (3.10) with the new definition for C and D together 
with the constraints ||(£ty),j | 2 = 1 , (t = l,...,n) . It gives n eigenvalues where each complex 
valued one appears together with its conjugate. With each real eigenvalue we can associate a real 
eigenfunction that solves a real eigenvalue problem like (2.2), while each complex one, together with 
its complex eigenvector, defines a pair of functions, corresponding to the real and the imaginary 
part of the complex eigenfunction associated with it, and a problem like (2.3). In an analogous 



7 


way we define a set of n real vectors, constructed from the eigenvectors of (3.10), which define an 
nxn real valued matrix whose rows are the above real n valued vectors and whose action on the set 
. . . ,<t> n ) defines a new set of functions which is a better approximation to the corresponding 
set of eigenvalue problems. 

4 Multigrid Algorithm 

From the above single grid method we construct a multigrid algorithm. The relaxation defined 
for the single grid algorithm involved a local step and a global step. The local step which uses 
the damped Jacobi relaxation has the property of damping errors that belong to non-dominant 
eigenfunctions on that grid. These functions are mainly high frequency functions which cannot be 
represented on a coarser grid. The local step serves therefore as a smoother for our problem. The 
global step updates the global quantities (the eigenvalues), as well as performing a projection of 
the solution onto an appropriate subspace that corresponds to the eigenvalues under consideration. 
This projection step is designed to control the important smooth eigenfunction to be calculated. 

In trying to construct a multigrid algorithm we have a few options in treating the global step. 
The first is to employ it on all levels of discretization. Such an approach is in contrast to a basic 
principle in multigrid; that is, global changes should be performed on coarse grids only [l]. Note 
that these global changes involve in most problems the convergence of the smooth components and 
therefore it makes sense to perform them on coarse grids. The effect of the global step on the 
non-smooth components is not important since these components should converge to zero as the 
algorithm proceeds. However, one has to make sure that the size of the non-smooth components is 
small relative to the relevant smooth ones before performing the global step, since small errors in 
the highly oscillatory components can cause large errors in the global quantities. This is another 
reason why it is preferred to perform the global step on coarse grids, where small errors in the high 
frequency functions have smaller effect on the outcome. 

4.1 Coarse Grid Correction 

For clarity assume first that M is spanned by a pair of complex conjugate eigenfunctions. In that 
case the relevant eigenvalue problem is a discretization of (2.3). The standard coarse grid equations 
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in a multigrid algorithm for that case using the FAS formulation are 


( A H - XI pi \ | 

(*») 

( w) ^ 

►-1 

1 

1 

{*" J 

„ W*) J 


II^IIh — r f > llV^II/r — 

where 

ThW 1 ) - A H I?w h - lf?A h w h 

r5 = rt + PM\%-\\^\\l 


( 4 . 1 ) 


( 4 . 2 ) 


( 4 . 3 ) 

( 4 . 4 ) 


where u > h is any fine grid function. 

The right hand side of the coarse grid equation depends heavily on the current approximation 
on the fine level, hence it would be different if a global step was carried out on the fine grid. In 
other words, the fine grid change 


4 > h 



would change the right hand side of the coarse grid problem to 


( 4 . 5 ) 


( ail <*12 \ ( T?{<t> h ) 

^ <*21 <*22 J \_ ^(xfi h ) 

Since the global step may change the solution on the fine grid significantly, and as a result also 
the right hand side of the coarse grid equations, it is better if the coarse grid right hand side 
takes this global change into account. The role of the equations H^ll# = r % , and ||V>*|||r = r yf 
is to prevent us from getting a zero solution to our problem. Any other equation that will have 
the same effect is equally good. Since the projection step leads to complicated equations for these 
constraints, we deviate here for these constraints, from the standard FAS. We basically ask that the 
coarse grid solution is of norm one for both the components under study. This leads to a significant 
simplification. This together with taking into account the global step as if done on the fine level 
leads us to the following coarse grid equations. 


( 4 . 6 ) 
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( A h - XI fil ^ ( an ai2 ^ f <t> H ^ _ ( an a i2 \ ( T?(<f> h ) \ ^ ^ 

^ -fil A h - XI J y a 2 1 a 22 J y tl> H J y a 2i a 22 J y T^(tp h ) J 

||an^ H + a n ip H \\ 2 H = 1, ||a 2 i^ ff + ot^ H \\]t ~ 1 ( 4 - 8 ) 

In this equation the unknowns are n and an, ai 2 ,a 2 i, a 22 . The iterative method for 

solving this problem will involve Jacobi relaxation for keeping the rest of the unknowns 

fixed, followed by a global step. The global step requires the solution of an eigenvalue problem like 
(3.10) with the following C and D 


C- ( < AH * H ~ > < AH ^ H ~ > 

y < A H <f> H - > < A H rp H - Tj?(i/> h ),ip H > 

D= ( <* tt A H > \ 

^ < tf> H , 1p H > < lf) H , 1p H > J 

A corresponding eigenvector defines the four a 1 a as in the single grid method. Having solved this 
problem approximately for <f> H , fi and an,ai 2 ,a 2 ia 22 the interpolation to be used is 



/ ^ 


W 

<t> k 


an a i2 J 


4 




^ a 2 i a 22 J ^ 


<t> k + {4> h - / A V) 

tj, h + ( ip H - ij?if> h ) 


\ 

J 


(4.11) 


We now come to the general case. 

dim 1/ = n. Assume that there are r complex eigenvalues and n-r real eigenvalues corresponding 
to this subspace. In this case the fine grid problem is of the form 


where £= (<£i, . . . ,<f>n) T , 


U k ~ L k )± h = 0 



V 


\ 


(4.12) 


(4.13) 

(4.14) 
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and 


{ Ai y.i 

-fii Ai 


\ 


Ar Mr 

-Hr A r 

Ar+1 


V 


An 


nxn 


Following the previous arguments, the coarse grid problem is then 


(4.15) 


U H - T)E£ h = EiZ(± h ) 


(4.16) 


IK^WI* = 1. (< =l,...,n) 


(4.17) 


where 


d(± h ) = 




(4.18) 


and E is an nxn real valued matrix, whose rows are obtained from a solution of the following 
eigenvalue problem 


det(C - pD) = 0, (4.19) 

where D is an nxn matrix whose elements are < 4f > > and C is an nxn matrix whose elements 

are < >• Since D is invertible there are n eigenvalues where complex ones 

appear in pairs of conjugates. A complex eigenvalue p correspond to a pair (A j,Hj) f°r some j < l 
while a real one corresponds to a real Xj for a j > l. The real and the imaginary part of a complex 
eigenvector will define two rows in E that correspond to the corresponding pair (A j,Hj)t while a 
real valued eigenvector will define one row in E in a place corresponding to that of the eigenvalue 
associated with it. 

Having solved the coarse grid problem approximately the following interpolation is to be used: 

t h - E{± h + ~ g4 h )}. (4.20) 
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4.2 Multigrid Algorithm 

Consider a sequence of grids < m) with mesh sizes h k satisfying 2 h k+1 = h k . Suppose on 
each grid operators A k are given in such a way that A k , (k < m) is an approximation to A k+1 . 
Assume also that interpolation operators I k _ v from coarse to fine grids , and a restriction operator 
I £ -1 , from fine to coarse grids, are given. The problem is to solve the following equation given on 
the finest level 


U»-A“)f* = o, ||/i(£“MI* = i (.' = 1 »)• (4.21) 

On levels (A: < m) the equations to be solved are of the form 

U* - A*)££* = ®ll(£* +1 ), l|/i(£*)(l| 2 = 1 (i=l n). (4.22) 

where 

lJk(fc£* +1 ) = £ k l!k+iU k+1 + I*+i{l^ 1 Gii fc+2 ) - £ m k+1 u k+1 } ( k < M ) (4.23) 

= 0 (4.24) 

where w fc+1 is any fine grid n-vector function. 

Note that the finest grid problem, i.e., for k=m, can be written as the general k-grid one with 
the choice E=I, the nxn identity matrix. 

Given an approximate solution {(^ m , rp m ) T , \,n, E} to the problem (4.22) the multigrid cycle 
for improving it is denoted by 


(**,**, A, /*,£) - MG{k, <f> k , rj> k , A, /i, E) (4.25) 

and is defined recursively as follows: 

if k=l then solve (4.22) by enough local relaxations together with global steps 
else 
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• Perform v\ local relaxation sweeps on (4.22), starting with (<f> k ,\p k , A, p,, E) and resulting 

in a new approximation (<^>*, t/> fc , A, p, E) 

• Starting with <f> k ~ l = ^ k ~ l = 7* -1 V> fc make 7 successive cycles of the type 

{4> k -\^ k -\\,H,E) «- MG(k — \,<j> k ~ l ,rl> k ~ l n,E) 

• Interpolate correction according to the formula <f> k = 4> k + i k k-M k ~ l - I ' c k~ 1 * k ) • 

• If k = m make the change 4> m «— E<t> m , then set E = I. 

• Perform 1^2 relaxation sweeps on (4.22) starting with {4> k , tp, n, E) and yielding 
( <f> k ,rp k X,n,E ), the final result of (4.25). 

FMG Algorithm 

In order to obtain full efficiency, the first approximation on a given level is obtained from a 
solution of the same problem on the next coarser level, which itself has been calculated in a similar 
way. The resulting algorithm is called (FMG) and is described next. Let n*_j be an interpolation 
operator (usually of higher order than Given the problem (4.21) , the N-FMG solution of 

that problem is: 

N-FMG Algorithm 
Set m = 1. 

Obtain initial approximation as in section (3.1). 

Calculate V’ 1 , A,/i) the solution of (4.21) for m=l by several relaxations, 
for m = 2 , . . . , M do: 

• Calculate <j> m <- t/> m <- II£_ 1 t/' m_1 

• Perform the cycle (^ m , rp m , A, p, E) <— MG(m, <£ m , rp m , A, p, E) N times. 


5 Numerical Examples 

Numerical experiments were carried out with the following operator 

_ ( A + c 2 a 0 ^ 
y bo A + c 2 j 
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in (0,l)x(0,l)> with ao x bo < 0. 

The eigenvalues of this operator are 


A ± tfi 

where A is an eigenvalue of A and ft = \/(-ao x to). 

Uniform discretization was used, where A was discretized using the standard 5-point Laplacian. 
n*+ J was bi-cubic interpolation, l£ +1 was bi-linear interpolation and l£ +1 was the 9-point full 
weighting operator. The following parameters were used in the multigrid cycling: = 2, = 

1,7 = 1. Initial approximation on the coarsest level was obtained as in section (3.1). Results for 
the 2-FMG algorithm are given in Tables 1-4. The L% norm of the residuals as well as the values 
of X h and pt* are given at the end of each cycle. Although a 2-FMG algorithm was used the results 
clearly show that 1-FMG algorithm solves the problem to the level of discretization errors. Observe 
that the imaginary part pt of the eigenvalue is exact on coarse grids since the truncation errors in 
approximating the off diagonal elements of A are zero, and only those elements contribute to the 
non symmetry of the problem. Also note the insensitivity of the results to the closeness of the 
eigenvalue to the imaginary axis. 
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level 

cycle 

||restdua/s ||2 

A* 

p h 

2 

1 

.279e+00 

-19.4809 

1.00000 


2 

.150e-01 

-19.4863 

1.00000 

3 

1 

.298e-01 

-19.6737 

1.00000 


2 

.151e-02 

-19.6757 

1.00000 

4 

1 

.312e-02 

-19.7228 

1.00000 


2 

.203e-03 

-19.7233 

1.00000 


Table 1: oq =1. bo = — 1. c 2 = 0. 


level 

cycle 

\\re 8 %duals \\2 

X h 

P h 

2 

i 

.279e+00 

-19.4809 

3.16227 


2 

.150e-01 

-19.4863 

3.16227 

3 

1 

.298e-01 

-19.6737 

3.16227 


2 

.151e-02 

-19.6757 

3.16227 

4 

1 

.312e-02 

-19.7228 

3.16227 


2 

.203e-03 

-19.7233 

3.16227 


Table 2: qq = 10. bo = — 1. c 2 = 0. 
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level 

cycle 

jjree»duols ||2 

X h 

BIS 

2 

i 

.279e+00 

0.51912 

3.16227 


2 

,150e-01 

0.51370 

3.16227 

3 

1 

.298e-01 

0.32634 

3.16227 


2 

.151e-02 

0.32421 

3.16227 

4 

1 

.312e-02 

0.27714 

3.16227 


2 

,203e-03 

0.27766 

3.16227 


Table 3: o 0 = 10. b 0 = -1. c 2 = 20. 


level 

cycle 

\\reeiduala \\2 

X h 


2 

i 

.279e+00 

0.51912 

1.00000 


2 

.150e-01 

0.51370 

1.00000 

3 

1 

.298e-01 

0.32634 

1.00000 


2 

.151e-02 

0.32421 

1.00000 

| 

1 

.312e-02 

0.27714 

1.00000 

■ 

2 

.203e-03 

0.27766 

1.00000 


Table 4: ao =1. 6q = — 1. c 2 = 20. 
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